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Abstract 

The computer algebra routines^ documented here empower you to 
reproduce and check many of the details described by an article on large 
deviations for slow-fast stochastic systems [6]. We consider a 'small' 
L| spatial domain with two coupled concentration fields, one governed by a 
. ^ 'slow' reaction-diffusion equation and one governed by a stochastic 'fast' 
linear equation. In the regime of a stochastic bifurcation, we derive two 
superslow models of the dynamics: the first is of the averaged model of 
the slow dynamics derived via large deviation principles; and the second 
is of the original fast-slow dynamics. Comparing the two superslow 
models validates the averaging in the large deviation principle in this 
parameter regime [6]. 
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1 Iterative computer algebra derives the model 

Construct a one element model of the 'slow' stochastic reaction diffusion 
equation 

ut = 9xxU + Asinu- (1 - 9xx)"^u - \/ea(1 - 9xx)"^ ct)(x, t) (1) 
such that u = at X = 0, 7t , 

near the deterministic bifurcation that occurs at A = 3/2 , to effects quadratic 
in the noise amplitude a. We seek the normal form where the evolution 
involves no convolutions [1, 5]. Also, transform the quadratic noise in the 
evolution. Throughout we adopt the Stratonovich interpretation of stochastic 
differential equations so that the ordinary rules of calculus apply. 

The stochastic slow model appears to be, when parameter A = j + A' and 
upon truncating the noise to just the first three sine modes, 

a = A'a - + lA')a3 + ^^a' - ^a(lct)i + i^a^4)3) + • • • (2) 

when the stochastic slow manifold is 

tl = d sin X + sin 3x 
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ea^^sinlxe lo * (j)2 + sin 3x e 5 -k (^7, 
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In outline, the algorithm iteratively determines the stochastic slow manifold 
model, then finally transforms to a weak model by replacing quadratic 
noises by their long time equivalents. Earlier research [2] explained the 
centre manifold rationale and the computational effectiveness of this simple 
algorithm, albeit there restricted to deterministic systems. 

» ssmaveq « 

7„ see cassmaveq.pdf for documentation 

« initialisation l>l> 

« linear noise effects » 

« quadratic noise effects \» 

sig: =small*sigma; 

let { small~6=>0 }; 

it:=l$ 

repeat begin 

« compute residual t>[> 

« update ssm » 

showtime ; 
end until res=0 or (it : =it+l) >20 ; 
write gssm:=sub(small=l,g) ; 
« transform quadratic noise » 
end; 



1.1 Initialisation 

Trivially improve printing. 

» initialisation « 
on div; off allfac; on revpri; 
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factor sigma, sin, small ; 
linelength 65$ 

Define the parameter A to be a small perturbation away from critical. Scale 
this with ordering parameter small in order to control truncation in the 
multiple small parameters. 

» initialisation «+ 
lamb:=3/2+small~2*lam; 

Linearise products of trigonometric functions via trigsimp. 

Define otm to be the decay rate of linear modes, here ara = — 3/2 + 
1 /(m^ + 1 ) , so that the spatial modes decay linearly like sin(mx)e^*"^^. 

I» initialisation «+ 
procedure alfa(m) ; (m"2-3/2+l/(m"2+l))$ 

Define the inverse of the linear operator, C^^ sin(mx) = sin(rax)/aTa , as the 
linear operator is £ = —3/2 — 9xx + (1 — 9xx)^^ • Note: we only define and 
use this for m > 2 . 

» initialisation «+ 

operator uinv; linear uinv; 

let uinv(sin(~m*x) ,xt)=>sin(m*x)/alf a(m) ; 

Define the linear operator, (1 — 9xx)^^ sin(rLx). 

» initialisation «1+ 

operator iddi; linear iddi; 
let { iddi(sin(~n*x) ,x) => sin(n*x) / (l+n"2) 
, iddi(sin(x) ,x) => sin(x)/2 }; 
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Paramterise solutions by an evolving amplitude a(t) (or 'order parameter'). 
Its evolution is dd/dt = d = g . 

» initialisation «+ 

depend a,t; 
let df(a,t)=>g; 

Then the most basic linear approximation to the dynamics on the element is 
u = dsinx where d = . Scale the amplitude to be small. 

» initialisation «+ 

u:=small*a*sin(x) ; 

v:=u/2; 

g:=0; 



1.2 Compute residual 

The parameter small conveniently controls the trunction in nonlinearity 
and other small parameters. The iteration terminates when the residual of 
the reaction diffusion equation is zero to the specified order of smallness. 
Note: this parametrisation with small should create deterministic models to 
errors 0(d^ + e^), or equivalent, as we scale e = eps with small~2. For some 
strange reason we need to do some operation on res in order for relevant 
terms to cancel — here I use trigsimp, but something else might serve. 

OO compute residual « 

sinu : =trigsimp (u-u"3/ 6+u"5/ 120-u" 7/5040 , combine) ; 
res : =-df (u,t)+df (u,x, 2)+lamb*sinu-iddi (u,x) 

-small*rooteps*sig*iddi (noise , x) ; 
res : =trigsimp(res , combine) ; 
write lengthres : =length(res) ; 
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Define ^/e which scales the size of the noise. Looks like we do not have to 
worry about rooteps/eps not simplifying. 

» initialisation «+ 
let rooteps~2=>eps ; 

1.3 Update the stochastic slow manifold 

Let T, tt, label the fast time of stochastic fluctuations so we can separate 
the stochastic fluctuations from the superslow evolution of the amplitude d. 
Also introduce xt to label both the subgrid spatial scales and time scales so 
we can group all factors in the space-time dynamics. 



depend tt , t ; 
depend x , xt ; 

Then update driven by the residual. Divide the evolution g by small to best 
keep track of the correct counting of the 'order' of a term. 

» update ssm « 

g : =g+ (gd : =secular (res , xt) ) /small ; 
u:=u+uinv(res-gd*sin(x) ,xt) ; 



1.4 Linear noise effects 

Introduce the noise in its spatial Fourier decomposition 



» initialisation «1+ 



oo 




n=1 
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Parametrise the amplitude with a. Truncate the spatial structure of noise. 
Three terms in the noise appears adequate to show the typical interactions 
between noise and other dynamics. Have not explored details of better 
resolution of the noise. 

» linear noise effects « 

operator phi; depend phi,tt,xt; 
noise:=for n:=l:3 sum phi(n,{})*sin(n*x) ; 

Let phi (n, {ml, . . .}) denote convolutions with exp(-ml*t) . . ., that is, 
45n,(m,,m2,...) = exp(-mi t) * exp(-Ta2t) -k ■ ■ ■ -k 4)n(t) ; 

so 

9t4>n,(m,,m2,..0 = 4'n,(m, ,m2,.-) + 4>n,(m2,.-) • 

But if we pull out a decay rate in 1 /e then keep bookkeeping correct by 
dividing by small"2 unless it is a sole convolution by ©(l/e). This latter 
case is only used in the next section. 

» linear noise effects «+ 

let { df (phi(~m,~p) ,t)=>df (phi(~m,~p) ,tt) 

, df (phi (~m, ~p) ,tt)=>( -first (p) *phi (m,p)+phi (m,rest (p) ) ) 

when deg(l/f irst (p) , eps)=0 
, df (phi (~m, ~p) ,tt)=>(-f irst(p)*phi(m,p)/small"2 

+phi (m, rest (p) ) /small" (if rest(p)=-[} then 1 else 2)) 

when deg(l/f irst(p) ,eps)=l 

>; 

Recall the equation for updates u' and g' is g' + C\i' = residual, where now 
the operator £ = 9t — 3/2 — 9xx + (1 — 9xx)^^ includes fast time variations. 
The operator secular extracts from the residual all those terms which would 
generate generate secular growth in the field u and so instead must be placed 
in the model's evolution g. The last rule here comes from integration by 
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parts and is essential in order to eliminate memory integrals (convolutions) 
in the model evolution. 

The if-clause in the last is needed for the next section to account for the 
convolution on the fast time only being of O (\/e) > ^^'^ so order increases 

» linear noise effects «+ 

operator secular; linear secular; 
let { secular (sin(~ni*x) ,xt)=>0 

, secular(sin(~ni*x)*~aa,xt)=>0 
, secular (sin(x) ,xt)=>l 
, secular(sin(x)*phi(~n,~p) ,xt)=> 
phi (n, {}) * (f or each r in p product (1/r)) 
*(if p neq{} and deg(l/f irst(p) ,eps)=l then small else 1' 

}; 



Extend the inverse operator to terms with fast time variations as well as 
fast (subgrid) space variations. Recursive procedure gungb extracts the non- 
secular parts of fluctuating sinx components. Have to adjust the smallness 
whenever the convolution transformed is on the e scale. 

» linear noise effects «+ 
procedure gungb (n,p) ; 
if p={} then else 
(gungb (n , rest (p) ) -phi (n , p) 

*(if deg(l/f irst(p) ,eps)=0 then 1 else small'2) 
) /first (p)$ 



let { uinv(sin(~m*x)*phi(~n,~p) ,xt)=>phi(n, (alfa(m)) .p)*sin(m*: 
, uinv(sin(x) *phi (~n, ~p) ,xt)=>gungb(n,p) *sin(x) 

}; 
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1.5 Quadratic noise effects 

Now let Zp denote multiple convolutions of in time of any term, zz(a,p) 
(though I only use Z for quadratic terms, it may well be able to replace the 
linear convolutions). That is, 

Z(m,,m2,...) = exp(-mit) *Z{TT^2,-) ^""^ 2() = 1 • 

» quadratic noise effects « 

operator zz; depend zz,tt,xt; 
let { zz(~a,{})=>a 

, df (zz(~a,~p) ,t)=>df (zz(a,p) ,tt) 

, df (zz(~a,~p) ,tt)=>-f irst (p) *zz(a,p)+zz(a,rest (p) ) 
when deg(l/f irst (p) , eps)=0 

>; 



To extract quadratic corrections to the evolution, use integration by parts 
so all non-integrable convolutions are reduced to the cannonical form of 
the convolution being entirely over one noise in a quadratic term, either 
<t>n<PmA...) or 4>n,(...)4>m- 

Have now made this very complicated for at least some of the cases when 
the convolutions may be over e-fast time scales. 

» quadratic noise effects «+ 

procedure gungd(n,p,m,q) ; 

if (p={»or(q={}) then phi(n,p)*phi(m,q) 

else if deg(l/f irst(p) ,eps)=deg(l/first(q) ,eps) then 

(gungd (n , rest (p) , m , q) +gungd (n , p , m , rest (q) ) ) 

/ (first (p) +f irst (q) ) 

*(if deg(l/f irst(p) ,eps)=0 then 1 else small"2) 
else if deg(l/f irst(p) ,eps)=l then 

(gungd(n,rest(p) ,m,q)*(if rest(p)={} then small else small' 
+gungd(n,p,m,rest (q) )*small"2 
) /first (p) *sub(rat=-f irst (q) /first (p) ,geom) 
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else gungd(m,q,n,p)$ 

let { secular(sin(x)*zz(~a,~p) ,xt) =>secular(sin(x)*a,xt) 
*(for each r in p product (1/r)) 
, secular(sin(~m*x)*zz(~a,~p) ,xt)=>0 

, secular (sin(x) *phi ("n, ~p) *phi (~m, ~q) ,xt) =>gungd(n,p,m,q! 
, secular(sin(x)*phi(~ii,~p)~2,xt) =>gungd(n,p ,n,p) 

}; 

Extend C^^ operator uinv to handle quadratic terms. First, integration by 
parts gives all integrable contributions from direct product terms. 

Have to similarly modify gunge for at least some of the cases when the 
convolutions may be over e-fast time scales. 

» quadratic noise effects «+ 

procedure gunge (n , p , m , q) ; 

if (p={»or(q={}) then 

else if deg(l/f irst(p) ,eps)=deg(l/f irst(q) ,eps) then 
( -phi (n , p) *phi (m , q) 

*(if deg(l/f irst(p) ,eps)=0 then 1 else sniall)''2 

+gunge(n,rest(p) ,m,q) 

+gunge(n,p,ni,rest(q)) 
)/ (first (p)+first(q)) 
else if deg(l/f irst(p) ,eps)=l then 
(-phi (n,p) *phi (ni,q) * small "2 

+gunge(n,rest(p) ,m,q)*(if rest(p)={]- then small else sma! 

+gunge(n,p,m,rest (q) )*small"2 

) /first (p) *sub(rat=-f irst (q) /first (p) ,geom) 
else gunge (m,q,n,p)$ 
let { uinv(sin(x) *phi (~n, ~p) *phi (~m, ~q) ,xt) 

=>gunge(n,p,m,q)*sin(x) 
, uinv(sin(x) *phi (~n, ~p) "2,xt)=>gunge(n,p,n,p)*sin(x) 

}; 
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Second, similar integration by parts gives integrable contribution from terms 
involving convolutions of products. 

» quadratic noise effects «+ 

procedure gungf(a,p); 
if p={} then else 

(gungf (a,rest(p))-zz(a,p))/first (p)$ 
let { uinv(sin(~l*x) *phi (~n, ~p) *phi (~m, ~q) ,xt) 

=>sin(l*x) *zz (phi (n,p) *phi (m, q) , {alf a(l) }) 
, uinv(sin(~l*x) *phi (~n, ~p) "2 ,xt) 

=>sin(l*x)*zz(phi(n,p) "2,{alfa(l)>) 
, uinv(sin(~l*x)*zz(~a,~p) ,xt)=>sin(l*x)*zz(a,alfa(l) .p) 
, uinv(sin(x)*zz(~a, ~p) ,xt)=>sin(x)*gungf (a,p) 
>; 



1.6 Transform quadratic noise 

Now proceed to transform the strong model to a weak model by replacing 
the quadratic noises by their effective long term drift and volatility. Earlier 
work [3] determines the rationale for the details of this transformation. 

Set small = 1 as it has done its job of truncating the nonlinear terms in the 
asymptotic expansion. 

» transform quadratic noise « 

small : =1 ; 

write "Now transforming the quadratic noises"; 

Now transform the quadratic noise into new noises y\) (psi) that are equivalent 
in their long time statistics: the operator long implements the long-time 
equivalent noises as determined earlier [3] . For now only transform up to two 
convolutions. These new noises have subscripts that uniquely identify them. 
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» transform quadratic noise «+ 

operator long; linear long; 
operator psi; depend psi,tt,xt; 
let { long(l,tt)=>l 

, long(phi(~i,{}) ,tt)=>phi(i,{}) 
, long(phi(~i,{})*phi(~j ,{~k» ,tt) 
=> 1/2* (if i=j then 1 else 0) 
+psi(i, j ,-[k»/sqrt(2*k) 
, long(phi(~i,{})*phi(~j ,{~k2,~kl» ,tt) 
=> (psi(i, j ,{kl})/sqrt(2*kl) 
+psi (i , j , {k2 ,kl}) /sqrt (2*k2) ) / (kl+k2) 

}; 

gg:=long(g,tt)$ 



Root sum squares of the determined noise coefficients; this procedure im- 
phcitly assumes that there is no correlation between the multitude of noises 
in these two terms in the amplitude equation. This assumption appears 
correct for these two terms in this PDE. In general we would need to do a 
QR factorisation of the noise terms. 

» transform quadratic noise «+ 

operator sumsqpsi; linear sumsqpsi; 

let { sumsqpsi (l,tt)=>0 

, sumsqpsi (psi (~i , ~j , ~p) ,tt)=>0 

, sumsqpsi (psi (~i , ~j , ~p) "2 ,tt)=>l 

, sumsqpsi (psi (~i , ~ j , ~p) *psi (~ii , ~ j j , ~pp) ,tt)=>0 

>; 



Have a look at the numerical coefficients. 

» transform quadratic noise «+ 
on rounded; print_precision 5; 

gg:=gg; 
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Extract the coefficients of the terms in a and a a, both mean and fluctuating. 



let abs(eps)=>eps; 

c20 : =sqrt (sumsqpsi (coef fn(coef fn(gg, sig, 2) ,a,0)"2,tt)) ; 
c21mean:=(coeffn(coeffn(gg,sig,2) ,a, 1) 

where psi (~i , ~ j , ~p)=>0) ; 
c21 : =sqrt (sumsqpsi (coef fn(coeffn(gg, sig, 2) ,a,l)~2,tt)) ; 

Switch back to the rational arithmetic mode for any other postprocessing. 

» transform quadratic noise «+ 
off rounded; 
showtime ; 

Executing the resultant code constructs the superslow model of the stochastic 
bifurcation in the slow averaged spdes. 

2 Model interacting fast-slow-superslow 
components 

This section constructs a superslow model of the fast-slow stochastic reaction 
diffusion equation 



near the deterministic bifurcation that occurs at A = 3/2 , to effects quadratic 
in the noise amplitude o", and seeks the normal form where the evolution 
involves no convolutions. Recall that throughout we adopt the Stratonovich 
interpretation of stochastic differential equations so that the ordinary rules 
of calculus apply. 



» transform quadratic noise «+ 



u-t = 9xxU- + A sin u — V 




(3) 
(4) 



such that u = v = Oatx = 0, 7t, 
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The presence of tlie three time scales in the dynamics — the fast v, the slow u, 
and the superslow evolution of the bifurcation amplitude a — means that 
the algorithm outlined here is one of the most technically challenging ones 
I have implemented. I recommend understanding simpler systems before 
attempting to understand the details of this section. 

The resulting model appears to be the following to some order in small 
parameters. In terms of the superslow evolving amplitude a(t), where 
u ~ asinx and v ~ j'lsinx , parameter A = I + A', and noise in just three 
sine modes, a stochastic differential equation for the amplitude is 



X'{]+ley)a-{l + l\' + le)a' + ^,a' 



16 ' 8" ' 64' 



+ 



6 4)36 5^*4)3 



6080 



(5) 



To errors 0[e), this evolution equation is identical to the slow model (2) of 
the averaged equation with fluctuations. 

The corresponding stochastic superslow manifold appears to be that the slow 
field 



u = asmx + 



Jga^ sin3x+ 1 



easmxe e 



*4)i 



^"v/eo" sin 2x 
^"v/ecrsin 3x 



27. 



e lo'*— e e 



_38t _10t 

e 5 -k — e e * 



4)2 

453 + 



whereas the fast field has 0(1) fluctuations 
V = 7 a sin X + jj^e sin 3x 



+ —= sin X 



2 2 2' 

[1 + le]e"e^*+ie"e^*e"e^* 



*1 
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a 

H — 1= sm 2x 



H — 1= sin 3x 
V e 



1 _I0+ 1 _38 



*2 

10, 



*3 



+ 



In outline, the algorithm iteratively determines the stochastic superslow 
manifold model [2, e.g.]. 

Seem to have to keep one or two orders higher in small as there is a division 
somewhere. So actually compute to residuals and errors one or two orders 
in small less than apparently allowed here — this fudges the computations 
so that enough is kept to do the cancellation, then later truncates. That is, 
small"6=>0 actually computes to errors ©(small'^). Here, for some reason 
we have to tread carefully to get up to fifth order terms in small: first, 
compute linear noise effects; then, second, when residuals are zero, up the 
order to retain quadratic terms in noise and continue iterating. Alternatively, 
we could just seek terms up to fourth order in small by let small~7=>0. 
But I do want to get to fifth order because of the challenge. 

» ssmuv « 
°/o see cassmaveq.pdf for documentation 
let { sigma"2=>0, small~8=>0 }; 
sig: =small*sigma; 
« initialisation >\> 
« linear noise effects l» 
« quadratic noise effects 00 
it :=1$ 

repeat begin 

<l<l update from fast residual l>t> 

if {resu,resv}={0,0} then clear sigma~2;yoimplicitly sigma''3=>0 
« update from slow residual » 
showtime ; 

end until {resu,resv}={0,0} and (sigma~2neq0) or (it:=it+l)>19 
"/owrite ussm:=sub(small=l,u) ; 
°/„write vssm:=sub(small=l,v) ; 



Tony Roberts, March 2, 2013 



2 Model interacting fast-slow-superslow components 



16 



Table 1: order of magnitude of convolution operators. 







1 








1 




1 


1 


1 


1 


e e ★ 


O 




0{e) 


0{e) 


0{e) 




o 


;e3/2) 










o 




0{e') 


0{e') 


0{e') 



write gssm:=sub(small=l,g) ; 

<l<l transform quadratic noise » 

end; 



2.1 Some initialisation things 

Define pra to be the relative decay rate of linear modes of the fast variable v 
on the element, here pra = m-'^ + 1 , so that the spatial modes in v decay 
linearly like sin(mx) exp(— pTuV^) • 

OO initialisation «+ 
procedure beta(m) ; (m"2+l)$ 

Define some of the inverse of a linear operator. Now there are significant 
subtleties here: each convolution with rates ©(l/e) are themselves. Thus 
smallness is hidden in the convolution rates; consequently we have to track 
them artificially though a parameter such as small. Use small to count 
both the direct e and the hidden ones in the convolutions, as well as the 
other small parameters. 

A further complication is that single bare convolution is actually 0[^/ej [4, 
equation (27)]; Table 1 lists the correct order of magnitude of various convo- 
lutions. This complication is simplified a little by separating convolutions 
that occur on the fast time from the convolutions that occur on the low time 
scale. 
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The linear equations for updates are 

— g sin X — Ut + Uxx + jU — v + Resu = , 
— evt + u + Vxx — V + Resv = . 

When considering mode sin mx, the Unear equations for updates are, including 
a correction g to the evolution only in the case of the critical ra = 1 , 

-g - Ut - (am - l/|3m)u-- v + Resu = 0, 
— evt + u — PraV + Resv = , 

for the previous defined constants cXm. and |3ra- The details of solving for 
updates are not critical to correctness of the results (as the results should only 
depend upon driving the residuals to zero), but the details will determine 
whether the iteration does converge to zero the residuals. 

Updates from the u equation For residuals of the slow-equation, make 
updates to the u-field driven by the u-residual, and correspondingly update 
the v-field in a way that will not change its v-residual at this order. We insist 
on not changing the v-residual because this update is considered second and 
we must not undo earlier corrections driven from the v-residual. Dividing the 
u-update by pra is sufficient for the v-update. List here first the deterministic 
updates, second the generic linear noise update, and last the updates for 
resonant terms. Use procedure gungb to extract the non-resonant parts 
of Resu. 

Account for smoothing effect of convolution on noise via the if-clauses. Table 1 
shows that when a term goes from multiple fast time convolutions to include 
one additional slow time convolution, then the order of the term increases 
by ^/e (the following provision assumes we linearise convolutions so that any 
one term only has convolutions on the same time scale). 

» linear noise effects «+ 

operator uuinv; linear uuinv; 
operator vuinv; linear vuinv; 
let { uuinv(sin(~ni*x) ,xt)=>sin(m*x)/alfa(ni) 
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, vuinv(sin(~ni*x) ,xt)=>sin(m*x)/alfa(m)/beta(m) 
, uuinv(sin(~m*x)*phi(~n,~p) ,xt) 
=> phi(n, (alf a(m) ) .p)*sin(m*x) 

*(if p neq{} and deg(l/first(p) ,eps)=l then small else 
, vuinv(sin(~m*x)*phi(~n,~p) ,xt) 

=> phi(n, (alf a(m) ) .p)*sin(m*x)/beta(m) 

*(if p neq{]- and deg(l/f irst(p) ,eps)=l then small else 
, uuinv(sin(x) *phi (~n, ~p) ,xt)=>gungb(n,p)*sin(x) 
, vuinv(sin(x) *phi (~n, ~p) ,xt)=>gungb(n,p)*sin(x)/beta(l) 

>; 



To deal with quadratic noise, the following appear to to be enough for errors 
no higher order than small^. That is, with the modifications made to gunge 
and gungd. We do not seem to need any extra transformations from the 
residual of the fast v-equation, probably because it is linear. 

» quadratic noise effects «+ 

let { uuinv(sin(~m*x) *phi (~n, ~p) *phi (~1 , ~q) ,xt) 

=> zz(phi(n,p)*phi(l,q) , {alf a(m) }) *sin(m*x) 

7o*(if p neq{} and deg(l/f irst(p) ,eps)=l then small elsi 

, vuinv(sin(~m*x) *phi (~n, ~p) *phi (~1 , ~q) ,xt) 

=> zz(phi(n,p)*phi(l,q) ,{alfa(m)})*sin(m*x)/beta(m) 
7o*(if p neq{} and deg(l/f irst(p) ,eps)=l then small elsi 

, uuinv(sin(~m*x)*zz(~n,~p) ,xt) => zz(n,alf a(m) .p)*sin(m*x! 

, vuinv(sin(~m*x)*zz(~n,~p) ,xt) 

=> zz(n,alfa(m) .p)*sin(m*x)/beta(m) 

, uuinv(sin(~m*x) *phi (~n, ~p) "2 ,xt) 

=> zz(phi(n,p) "2 , {alf a(m) }) *sin(m*x) 

°/o*(if p neq{]- and deg(l/f irst(p) ,eps)=l then small elsi 

, vuinv(sin(~m*x) *phi (~n, ~p) ~2 ,xt) 

=> zz(phi(n,p) "2 , {alf a(m) }) *sin(m*x) /beta(m) 

7o*(if p neq{} and deg(l/f irst(p) ,eps)=l then small elsi 

, uuinv(sin(x) *phi (~n, ~p) *phi (~1 , ~q) ,xt)=>gunge(n,p,l,q)*s: 

, vuinv(sin(x) *phi (~n, ~p) *phi (~1 , ~q) ,xt) 
=> gunge(n,p,l,q)*sin(x)/beta(l) 
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, uuinv(sin(x) *phi (~ii, ~p) "2 ,xt)=>gunge (n,p,n,p) *sin(x) 

, vuinv(sin(x) *phi (~n, ~p) "2 ,xt)=>gunge (n,p,n,p) *sin(x) /bet; 

, uuinv(sin(x) *zz(~n, ~p) ,xt)=> 

( uuinv(sin(x)*zz(n,rest(p)) ,xt) 
-zz(n,p) *sin(x) )/first(p) 
, vuinv(sin(x) *zz(~n, ~p) ,xt)=> 

( vuinv(sin(x)*zz(n,rest(p)) ,xt) 
-zz(n,p) *sin(x) ) /first (p) /beta(l) 

}; 



Updates from the v equation Corrections to the u and v fields arise 
from the v-residual. However, because we consider this residual first in each 
iteration (not that first makes a lot of sense in an iterative loop), we are 
free to modify field u in a way that would affect the residuals at the same 
order. The key aspect is that we must not affect the residuals at a lower 
order in the u-residual — achieving this aspect is hard enough, which is why 
I implement corrections from the v-residual first. 

First define the deterministic updates. 

» linear noise effects «+ 

operator vvinv; linear vvinv; 
operator uvinv; linear uvinv; 

let { vvinv(sin(~ni*x) ,xt)=>sin(m*x)*(alfa(m)-l/beta(m)) 
/alf a(m) /beta(m) 
, uvinv(sin(~m*x) ,xt)=>-sin(m*x) /alf a(m) /beta(m) 
, vvinv(sin(x) ,xt)=>sin(x)/beta(l) 
, uvinv(sin(x) ,xt)=>0 



Second deal with the variety of linear noise terms. When a noise term in the 
residual is a convolution over the slow-scale, then the convolution is smooth 
and its time derivative correspondingly of the same order so that the evt 
causes no problem. 
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» linear noise effects «+ 

, vvinv(sin(~m*x)*phi(~n,~p) ,xt) 
=> phi(n,p)*sin(m*x)/beta(m) 
when p neq {} and deg(l/f irst(p) ,eps)=0 

, uvinv(sin(~m*x)*phi(~n,~p) ,xt) => 

when p neq {} and deg(l/f irst(p) ,eps)=0 

The critical mode is no different when forcing in the v-residual. 

» linear noise effects «+ 

, vvinv(sin(x) *phi (~n, ~p) ,xt) 
=> phi (n,p) *sin(x) /beta(l) 
when p neq {} and deg(l/f irst(p) ,eps)=0 

, uvinv(sin(x) *phi (~n, ~p) ,xt) => 

when p neq {} and deg(l/f irst(p) ,eps)=0 

However, when the noise term in the v-residual is not smooth, either because 
it is a bare white noise or because it is a convolution over the fast time scale, 
then we must be more careful because the evt term is important. Because 
the update has to be relatively large, ^ we have to use the u-field to cancel 
the effect in the u-residual of updates from the v equation. 

Here if the convolution is the first fast-time convolution, p={}, then choose 
the correct scale in small as then the convolution is only O^i/e). 

» linear noise effects «+ 

, vvinv(sin(~m*x)*phi(~n,~p) ,xt) 

=> phi(n, (beta(m)/eps) .p)*sin(in*x)/eps 

/ (if p={} then small else 1) 

when p={} or deg(l/f irst(p) ,eps)=l 

, uvinv(sin(~m*x)*phi(~n,~p) ,xt) 

=> phi(n, (beta(m)/eps) .p)*sin(m*x)/beta(m) 
*(if p={} then small else small~2) 

conjecture that it is this that affects the management of the smallness parameter. 
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when p={} or deg(l/f irst(p) ,eps)=l 



The critical mode is no different when forcing in the v-residual. 

» linear noise effects «+ 

, vvinv(sin(x) *phi (~n, ~p) ,xt) 

=> phi(n, (beta(l) /eps) .p)*sin(x)/eps 

/ (if p={} then small else 1) 

when p={} or deg(l/f irst(p) ,eps)=l 

, uvinv(sin(x) *phi (~n, ~p) ,xt) 

=> phi(n, (beta(l) /eps) .p)*sin(x)/beta(l) 
*(if p={} then small else small"2) 
when p={} or deg(l/f irst(p) ,eps)=l 



Second deal with a variety of quadratic noise terms. When a noise term 
in the residual is a convolution over the slow-scale, then the convolution is 
smooth and its time derivative correspondingly of the same order so that the 
evt causes no problem. 

» quadratic noise effects «+ 

let { vvinv(sin(~m*x)*zz(~n,~p) ,xt) 

=> zz(n,p)*sin(m*x)/beta(m) when deg(l/f irst(p) ,eps)=0 
, uvinv(sin(~m*x)*zz(~n,~p) ,xt) => when deg(l/f irst (p) , e] 
, vvinv(sin(~m*x) *phi (~n, ~p) "2 ,xt) 

=> phi (n,p) "2*sin(m*x) /beta(m) 

when p neq{} and deg(l/first(p) ,eps)=0 
, uvinv(sin(~m*x) *phi (~n, ~p) ~2 ,xt) => 

when p neq{]- and deg(l/f irst(p) ,eps)=0 
, vvinv(sin(~m*x) *phi (~n, ~p) *phi (~1 , ~q) ,xt) 

=> phi (n,p) *phi (1 , q) *sin(m*x) /beta(m) 

when p neq{]- and q neq{} 

and deg(l/f irst (p) ,eps)+deg(l/f irst(q) ,eps)=0 
, uvinv(sin(~m*x) *phi (~n, ~p) *phi (~1 , ~q) ,xt) => 
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when p neq{} and q neq{} 

and deg(l/f irst (p) ,eps)+deg(l/f irst(q) ,eps)=0 



The critical mode is no different when forcing in the v-residual. 

» quadratic noise effects «+ 

, vvinv(sin(x) *zz(~n, ~p) ,xt) 

=> zz(n,p)*sin(x)/beta(l) when deg(l/f irst(p) ,eps)=0 
, uvinv(sin(x) *zz(~n, ~p) ,xt) => when deg(l/f irst(p) ,eps)= 
, vvinv(sin(x) *phi (~n, ~p) "2 ,xt) 

=> phi(n,p) '2*sin(x)/beta(l) 

when p neq{} and deg(l/first(p) ,eps)=0 
, uvinv(sin(x) *phi (~n, ~p) "2 ,xt) => 

when p neq{]- and deg(l/f irst(p) ,eps)=0 
, vvinv(sin(x) *phi (~n, ~p) *phi (~1 , ~q) ,xt) 

=> phi (n,p) *phi (1 ,q) *sin(x) /beta(l) 

when p neq{} and q neq{} 

and deg(l/f irst (p) ,eps)+deg(l/f irst(q) ,eps)=0 
, uvinv(sin(x)*phi(~n,~p)*phi(~l,~q) ,xt) => 
when p neq{} and q neq{} 

and deg(l/f irst (p) ,eps)+deg(l/f irst(q) ,eps)=0 

}; 



Critical: linearise convolutions over different time scales I contend 
that we also want to simplify the convolutions because convolutions of the fast 
time scale e are qualitatively different from convolutions over slow time scales. 
Thus we do the following 'linearisation' of convolutions: whenever the first 
two convolutions are over different time scales, we transform the convolution 
into the sum of two convolutions. Change of variables in integration shows 
that 



1 



a 



-at 



-(3t^ 
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I have not used this transform in other apphcations because of the necessity 
to avoid division by zero when a = (3 Here, we are concerned with 
convolutions over different time scales and so use this formula where, for 
example, rate |3 is replaced by fast rate |3/e: 



-at -J^t 1 
e * e e ★ 



e/|3 



1 - ea/(3 



Thus we need to divide by 1 — r for various r = ea/(3 so store its power 
series in the variable geom, with small to account for the powers of e. The 
correctness of the following transformation is critical. 

» linear noise effects «+ 

geoin:=for n:=0:deg((l+small"2)"9,small)/2 sum (rat*small"2) "n$ 
let { phi(~n,~p) => (phi(n,f irst(p) .rest(rest(p))) 

-phi (n, second (p) . rest (rest (p) ))* (if rest (rest (p))={} 
or degd/first (rest (rest (p) ) ) ,eps)=l 
then small else small"2) 
) *sub(rat=f irst (p) /second (p) ,geom) /second (p) 
when length(p)>l and deg(l/f irst (p) ,eps)=0 
and deg(l/second(p) ,eps)=l 
, phi(~n,~p) => (phi(n,second(p) .rest(rest(p))) 

-phi (n, first (p) . rest (rest (p) ) ) * (if rest (rest (p) )={} 
or degd/first (rest (rest (p) ) ) ,eps)=l 
then small else small~2) 
) *sub(rat=second(p) /first (p) ,geom) /first (p) 
when length(p)>l and deg(l/f irst (p) ,eps)=l 
and deg(l/second(p) ,eps)=0 

>; 



^However, in general, maybe I should do this linearisation in order to reduce expressions 
to a more canonical form. 
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2.2 Update from residuals of the fast equation 

The parameter small, controls the truncation in nonhnearity and in smaU 
parameters. The iteration terminates when the residual of the reaction 
diffusion equation is zero to the specified order of nonlinearity. For some 
reason we need to do something nontrivial to the residual in order to force 
cancellation of terms so I apply trigsimp. Note the multiplication and 
division by small in order to cater for other divisions by small affecting the 
error truncation. 

» update from fast residual « 

resv: =-small''2*eps*df (v,t)+df (v,x, 2) -v+u 

+small*rooteps*sig*noise ; 
resv: =trigsimp(small~2*resv, combine) /small~2; 
write lengthresv:=length(resv) ; 
u:=u+uvinv(resv,xt) ; 
v:=v+vvinv(resv,xt) ; 



2.3 Update from residuals of the slow equation 

Similarly update from the residual of the slow equation. Divide the evolution g 
by small to best keep track of the correct counting of the 'order' of a term. 

» update from slow residual « 

sinu : =trigsimp (u-u"3/ 6+u"5/ 120-u" 7/5040 , combine) ; 
resu: =-df (u, t)+df (u,x, 2)+lamb*sinu-v; 
resu: =trigsimp(small~2*resu, combine) /small~2 ; 
write lengthresu:=length(resu) ; 

g : =g+ (gd: =secular (small "2*resu,xt) / small "2) /small ; 
u: =u+uuinv(resu-gd*sin(x) ,xt) ; 
v:=v+vuinv(resu-gd*sin(x) ,xt) ; 
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Executing the resultant code constructs the superslow model of the stochastic 
bifurcation in the fast-slow system of spdes. 
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